/*
    Copyright 2006-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/** \file 

    Contains declaration of the molecule_grain_model class. */

// $Id$

#ifndef __molecule_grain_model__
#define __molecule_grain_model__

#include "grain.h"
#include "thermal_equilibrium_grain.h"
#include "boost/shared_ptr.hpp"
#include "preferences.h"
#include "grain_model.h"

namespace mcrx {
  template <template<class> class, typename> class molecule_grain_model;
};


/** This class is a grain model consisting of a molecule of certain
    density and opacity. For now, it just assumes that the opacity is
    independent of temperature, so at the moment there's really no
    difference between this and the grain_collection_grain_model. */
template <template<class> class chromatic_policy, 
	  typename T_rng_policy = mcrx::global_random> 
class mcrx::molecule_grain_model : 
  public mcrx::grain_model<chromatic_policy, T_rng_policy> {
protected:
  /// The grain object that holds the cross section data.
  boost::shared_ptr<emitting_grain> grain_;
  /// The mass of the molecule, in the units of the grain object
  T_float mass_; 

  /** Calculates the extinction curve of the grain model. */
  void calculate_extinction_curve();

public:
  molecule_grain_model(const string& filename, const Preferences& p, 
		       const T_unit_map& units);
  molecule_grain_model (const molecule_grain_model& g) :
    grain_model<chromatic_policy, T_rng_policy>(g),
    grain_(g.grain_->clone()),
    mass_(g.mass_) {};
  virtual boost::shared_ptr<grain_model<chromatic_policy, T_rng_policy> > 
  clone() const {
    return boost::shared_ptr<grain_model<chromatic_policy, T_rng_policy> >
      (new molecule_grain_model(*this));};

  const array_1& alambda() const {return grain_->alambda();};
  const array_1& elambda() const {return grain_->elambda();};
  /** Resamples the grain object to use the specified wavelength
      vector. Also updated the extinction curve so the two are
      consistent. Note that this function is NOT necessarily safe to
      call several times with the same wavelength vector, because
      resampling with identical points does not necessarily preserve
      the data. */
  void resample (const array_1& lambda, const array_1& elambda=array_1()) {
    grain_->resample(lambda,elambda);
    calculate_extinction_curve();
  };

  /** This redefinition of the virtual scatterer::set_wavelength
      function just calls the resample function (which also
      recalculates the extinction curve). */
  virtual void set_wavelength (const array_1& lambda) {
    resample(lambda);
  };

  /** This gets called by the template virtual hack in
      grain_model::calculate_SED. */
  template<typename T>
  void calculate_SED_virtual(const blitz::ETBase<T>& intensity, 
			     const array_1& m_dust, array_2& sed,
			     array_2* temp) const;

  /** This function returns the total absorption cross section for the
      dust model given a certain dust mass. Intended for testing. */
  array_1 sigma_tot(T_float m_dust) const {
    using namespace blitz;
    array_1 st(grain_->sigma().extent(secondDim));
    st= grain_->sigma()(0, tensor::i)*m_dust;
    return st;
  };

};

template <template<class> class chromatic_policy, 
	  typename T_rng_policy> 
mcrx::molecule_grain_model<chromatic_policy, T_rng_policy>::
molecule_grain_model(const string& filename, const Preferences& p, 
		     const T_unit_map& units) :
  grain_model<chromatic_policy, T_rng_policy>()
{
  const std::string data_directory =
      word_expand(p.getValue("grain_data_directory", std::string ()))[0];
  const std::string full_path = data_directory + "/" + filename;

  // copy units from those supplied into *scatterer* unit map.
  T_unit_map& u = 
    mcrx::scatterer<chromatic_policy, T_rng_policy>::units_;
  u["length"] = units.get("length");
  u["wavelength"] = units.get("wavelength");
  u["mass"] = units.get("mass");

  // The grains use their preferred units, which are also copied into
  // the grain model units.
  grain_.reset(new thermal_equilibrium_grain(full_path, p));

  T_unit_map& u2 = 
    mcrx::grain_model<chromatic_policy, T_rng_policy>::units_;
  u2["length"]=grain_->units().get("length");
  u2["wavelength"]=grain_->units().get("wavelength");
  u2["mass"] = units.get("mass");

  // now read the molecular mass from the file
  CCfits::FITS f(full_path, CCfits::Read);
  open_HDU(f,"AXES").readKey("MOLMASS",mass_);
}
  

template <template<class> class chromatic_policy, 
	  typename T_rng_policy> 
void mcrx::molecule_grain_model<chromatic_policy, T_rng_policy>::
calculate_extinction_curve()
{
  using namespace blitz;
  using namespace boost;

  std::cout << "Calculating extinction curve" << std::endl;

  // we need to convert units. The extinction curve should be in the
  // units given by the scatterer units object, normalized to 1 mass
  // unit in that system. In the calculations below, the grain
  // quantities need to be converted. sigma has units of length^2.
  // We also need to convert our mass units to the scatterer mass units.
  const T_unit_map& us = 
    mcrx::scatterer<chromatic_policy, T_rng_policy>::units();
  const T_float length_conversion = 
    units::convert(grain_->units().get("length"), 
		   us.get("length"));
  const T_float mass_conversion = 
    units::convert(this->units().get("mass"), 
		   us.get("mass"));

  // The units of the extinction curve is area/mass, ie we should
  // divide by the mass conversion.

  // we don't currently convert the wavelength vector, though we want
  // to check
  assert(grain_->units().get("wavelength")==
	 us.get("wavelength"));

  // calculate extinction curve. this is based on the *emission*
  // wavelength grid also, since we know that there is only one
  // "species" in the opacity data, we don't have to sum
  const int n = grain_->elambda().size();
  array_1 abs(n);
  abs = grain_->esigma() (0, Range::all());
  abs *= length_conversion*length_conversion/(mass_*mass_conversion);

  // scattering opacity is assumed to be zero
  array_1 sca(n);
  sca = 0;

  array_1 g(n);
  g = 0;

  ASSERT_ALL(abs>=0);

  // CHECK THAT THESE VALUES ARE OK

  // now copy these to the scatterer members.
  assign(this->opacity_, (abs+sca));
  assign(this->albedo_, sca/(abs+sca));
  this->phfunc.set_g (g);

  DEBUG(1,								\
	ofstream crap("extinction_curve.txt");				\
	for(int i=0;i<grain_->alambda().size();++i)			\
	  crap << grain_->alambda()(i) << '\t' << (abs(i)+sca(i)) << '\t' \
	       << sca (i)/(abs (i) +sca (i)) << '\t'			\
	       << g (i) << endl;					\
	);
}


template <template<class> class chromatic_policy, 
	  typename T_rng_policy> 
template <typename T>
void
mcrx::molecule_grain_model<chromatic_policy, T_rng_policy>::
calculate_SED_virtual(const blitz::ETBase<T>& intensity, 
		      const array_1& mass,
		      array_2& sed, array_2* temp) const
{
  using namespace blitz;

  if(temp)
    // temp does not make sense for this grain
    *temp = blitz::quiet_NaN(T_float());
    
  // dn is the number of molecules per unit mass, so 1/mass_
  array_1 dn(1);
  dn=1./mass_;

  grain_->calculate_SED_from_intensity(intensity,
				       mass,
				       dn,
				       sed,
				       true, 0);
}


#endif
